Langevin agglomeration of nanoparticles interacting via a central potential 
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Nanoparticle agglomeration in a quiescent fluid is simulated by solving the Langevin equations of 
motion of a set of interacting monomers in the continuum regime. Monomers interact via a radial, 
rapidly decaying intermonomer potential. The morphology of generated clusters is analyzed through 
their fractal dimension df and the cluster coordination number. The time evolution of the cluster 
fractal dimension is linked to the dynamics of two populations, small (k < 15) and large (k > 15) 
clusters. At early times monomer-cluster agglomeration is the dominant agglomeration mechanism 
(df — 2.25), whereas at late times cluster-cluster agglomeration dominates (df = 1.56). Clusters are 
found to be compact (mean coordination number ~ 5), tubular, and elongated. The local, compact 
structure of the aggregates is attributed to the isotropy of the interaction potential, which allows 
rearrangement of bonded monomers, whereas the large-scale tubular structure is attributed to its 
relatively short attractive range. The cluster translational diffusion coefficient is determined to be 
inversely proportional to the cluster mass and the (per-unit-mass) friction coefficient of an isolated 
monomer, a consequence of the neglect of monomer shielding in a cluster. Clusters generated by 
unshielded Langevin equations are referred to as ideal clusters because the surface area accessible to 
the underlying fluid is found to be the sum of the accessible surface areas of the isolated monomers. 
Similarly, ideal clusters do not have, on average, a preferential orientation. The decrease of the 
numbers of clusters with time and a few collision kernel elements are evaluated and compared to 
analytical expressions. 

PACS numbers: 61.43.Hv, 82.70.-y, 47.57.cb, 47.57.J- 



I. INTRODUCTION 

Nanoparticle aggregates are of significant importance in technological and industrial processes such as, for example, 
combustion, filtration, and gas-phase particle synthesis. In addition, colloidal aggregates play an important role in, 
for example, the pharmaceutical industry, food processing, paintings, and polymers. The fractal nature of these 
aggregates has profound implications on their transport [1] and thermal [2] properties. Fractal aggregates arise from 
the agglomeration of smaller units, herein taken to be spherical and referred to as monomers, that do not coalesce but 
rather retain their identity in the resulting aggregate. In a quiescent fluid the main mechanism driving agglomeration 
is diffusion: accordingly individual monomers may be modelled as interacting Brownian particles whose motion 
and dynamics is described by a set of Langevin equations. Langevin simulations have been used to study aggregate 
formation [3] , aggregate collisional properties [4] , the limits of validity of the Smoluchowski equation [5] , and aggregate 
films [6]. 

In this study we investigate nanoparticle agglomeration and the diffusive motion and growth of the resulting 
aggregates by relying solely on the Langevin equations of motion of a set of interacting monomers in three dimensions. 
The monomer-monomer interaction potential is taken to be either a model potential, composed of a repulsive and a 
short-ranged attractive part, or a Lennard- Jones intermolecular potential integrated over the monomer volumes. Both 
potentials are spherically symmetric and rapidly decaying. Spherical symmetry implies that intermonomer forces are 
central: no angular force limits reorientation of an attached monomer. We shall argue that these two features have 
profound implications on the small and large-scale structure of the generated aggregates. Unlike the works mentioned 
above, no assumptions are made about the structure (e.g., by specifying the cluster fractal dimension) or the mobility 
of the aggregates. Our approach is reminiscent of other applications of Langevin simulations in dilute colloidal 
suspensions to study diffusion-induced agglomeration and the structures it gives rise to [7] . Stochastic agglomeration 
(particle motion and deposition) has also been discussed in terms of Brownian dynamics, see, for example, Refs. [8, 9]. 

An inherent difficulty of the adopted mesoscopic description of nanoparticle agglomeration, a description whereby 
the effect of collisions of fluid molecules with larger solid nanoparticles is modelled by a random force, is the choice 
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of the stochastic properties of the random force. These properties, in particular the noise strength, are usually 
specified by invoking the Fluctuation Dissipation Theorem (FDT). We will use the FDT to relate the amplitude of 
the fluctuations of the random force acting on a monomer to its friction coefficient. In applying the FDT, however, 
monomer shielding will be neglected in that the friction coefficient of a monomer in an aggregate will be taken to be 
independent of its state of aggregation. This approximation of the hydrodynamic forces acting on a monomer has 
been referred to as the free draining approximation [10]. The monomer friction and diffusion coefficients are related 
to the monomer surface area accessible (or exposed) to surrounding fluid molecules [11]. The accessible surface area is 
the fraction of the geometric surface area that is active in momentum and energy transfer from the underlying fluid to 
the monomer: it, thus, determines the monomer and aggregate transport properties. We will argue that the accessible 
surface area of a generated fc-monomer cluster (also referred to as a fc-mer), is k times the accessible surface area of 
an isolated monomer, i.e., the accessible surface area of a sphere. We shall, thus, refer to the clusters generated by 
unshielded Langevin equations as ideal clusters with respect to their dynamic properties. Monomer shielding may, 
equivalently, be described in terms of the cluster shielding factor defined as the ratio of the average (per-unit-mass) 
friction coefficient of a monomer in an aggregate to the (per-unit-mass) friction coefficient of an isolated monomer. 

In a non-ideal cluster the monomer accessible surface area decreases due to shielding leading to a decrease of the clus- 
ter friction coefficient, and a consequent increase of its diffusion coefficient. Filippov et al. [2] presented an expression 
for the shielding factor of a cluster by comparing the total heat transfer to an aggregate to the product of the number 
of monomers in the aggregate times the heat transfer to an isolated monomer (under the same conditions) . Monomer 
shielding may be understood by noting that the fluid concentration boundary layers of neighbouring monomers in a 
cluster overlap non-additively, thereby the fluid molecule-monomer collision rate decreases. In a Langevin description 
of cluster diffusion the change of the monomer accessible surface, and hence the change of the cluster shielding fac- 
tor, renders the strength of the stochastic thermal noise a time-dependent function of the local arrangement of each 
monomer. We do not attempt such a modification of the noise strength, adopting the form of the FDT that leads to 
unshielded Langevin equations. 

Since we solve numerically the Langevin equations for interacting monomers rather than the aggregate equations of 
motion, information on the dynamics of aggregate formation has to be inferred from the simulation output. We obtain 
this datum, together with the detailed structure of each aggregate, from the record of the collisions it underwent and 
its eventual restructuring, using techniques borrowed from graph theory [12]. 

The paper is organized as follows: Section II provides the theoretical framework for the Langevin simulations. 
Emphasis is placed on the description and justification of the monomer-monomer potentials used in the numerical 
experiments. Section III describes the numerical method, and it introduces the quantities monitored to investigate the 
statics and dynamics of the generated aggregates. Section IV presents the static properties of the aggregates, namely 
their fractal dimension, their morphology (specifically, the cluster coordination number), and some indications of 
cluster restructuring. Section V summarizes the cluster dynamic properties, in particular, their translational diffusion 
coefficient, their mean Euler angles, the decay of the total number of clusters with time, and the determination of 
a limited number of agglomeration kernel elements. The final remarks in Sec. VI summarize the main results and 
conclude the paper. 



II. MODEL DESCRIPTION 
A. Monomer Langevin equations of motion 

We investigate the non-equilibrium dynamics of interacting clusters in the continuum regime (fluid mean free path 
smaller than the monomer diameter) by solving the Langevin equations of motion of a dilute system of N interacting 
monomers in a quiescent fluid in three dimensions. We use the word cluster with the same meaning as the term 
aggregate in Ref. [13], i.e., a set of physically bound spherules (monomers). The i-th monomer obeys the Langevin 
equation 

mir i = F i -/3imir < +W < (t) , i = l,...N , (1) 

where mi in the monomer mass, its position in three-dimensional space, Fj an external force, /3i the monomer 
friction coefficient per unit mass, and W a noise term that models the effect of collisions between the monomer and 
the molecules of the surrounding quiescent fluid. We consider that each monomer feels a Stokes drag (continuum 
regime), the friction coefficient per unit mass being j3\ = 37r/i/<7 /mi where (if is the fluid viscosity, and a the monomer 
diameter. Henceforth all friction coefficients will be defined per unit mass. The monomer friction coefficient may also 
be expressed as the inverse monomer relaxation time (3\ = Tj~ where T\ = p p cr 2 / (18 (if) with p p the monomer material 
density. As argued in the Introduction, the use of Stokes drag implies that the monomer surface area accessible to fluid 
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molecules equals the accessible surface area of a monomer irrespective of its state of aggregation. The implications 
of this approximations are explored in Sec. VA. The noise is assumed to be Gaussian and white (delta-correlated in 
time) 

(W?(t))=0 , and (W l k (t)W l j (t'))=T8 ij 6 kl S(t-t') , (2) 

where angular brackets (. . .) denote an ensemble average over realizations of the random force, subscripts denote 
monomer number (i,j — l,...N), and superscripts Cartesian coordinates (k,l = x,y,z). The noise strength T is 
determined from the Fluctuation Dissipation Theorem applied to each monomer: it evaluates to T = 2/3imifcsT with 
Ub Boltzmann's constant and T the system temperature [14]. 

The force acting on the i-th monomer arises from its interaction with all the other monomers. It is considered to 
be conservative 

F i = -V P( I7 j , (3) 

where Ui is the total intermonomer potential the i-th monomer feels. We assume that it derives from a pair-wise 
additive, two-body, radial potential Uijfrij) 

N 

where the radial distance is = |r^ — r j | - 

Equations (1) may be cast in dimcnsionless form, a form more convenient for their numerical solution. By choosing 
length I*, time t* , and mass m* characteristic scales we introduce the dimensionless variables, denoted by tilde, 

r = l*r , t = T*t , mi = m*rhi . (5) 

These characteristic scales fix the characteristic system temperature T* to 



Equations (1) and (2) in dimcnsionless and component-wise form become 
d 2 f' 1 dU - df\ 1 T ~ 

i * _|_ 

dt 2 mi dr\ di fh\ 



^ - + —W!(i) ; l = x,y,z ; i = l,...N , (7a) 



(W7(t))=0 , and (W!®W?(?))=2 1 m 1 T5 ij 6 k i5(t-F) , (7b) 
where we introduced the dimensionless variables 

A -At* , f-£ , , Wl^Wl . (8) 

Equations (7) show that three independent dimcnsionless variables (mi, $i, T) determine the dynamics of the system. 
A natural, but not unique, choice of units is m* = mi, r* = n, and I* = a, where a is a characteristic length scale of 
the intermonomer potential (to be taken to be the monomer diameter). This choice leads to Pi = rhi = 1 in Eq. (7a), 
a form of the Langevin equations we will use in our numerical simulations. All simulations were performed at T = 0.5. 

The characteristic scales may be evaluated for a typical case of agglomeration of combustion-generated nanoparticles. 
A typical soot monomer of material density p p — 1.3g/cm 3 and characteristic size I* = a = 20nm has a characteristic 
temperature T* 

18 2 7TU 2 f C7 

T* = 7*7 — 650K , (9) 
6k B p P 

when suspended in air at 300K of dynamic viscosity /U/ = 1.85 • 10 _4 g/(cm-sec). Therefore, T = 0.5 corresponds to 
approximately T ~ 300K. The Stokes monomer relaxation time, the characteristic time scale, is t\ ~ 1.6 x 10 _9 sec. 
Heat, mass, and momentum transfer between particles and the carrier gas depend on the Knudsen number Kn = 
2X g /a, where X g is the carrier gas mean free path: for Kn <C 1 these transfer processes occur in the continuum regime. 
The air mean free path at atmospheric pressure and 293K is A a i r = 66nm. Hence, our simulations are appropriate 
either for aerosol agglomeration at high pressures (A g <~ p~ with p g the carrier gas pressure) or for agglomeration 
of non-charged colloids in liquids. 
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B. Monomer-monomer interaction potential 



The monomer-monomer potential is chosen to mimic interaction of hard-core monomers sticking upon collision. As 
such it will be taken to be rapidly decaying. The effect of the range of repulsive interactions in two-dimensional colloidal 
aggregation was extensively discussed in Ref. [16], whereas Ref. [7] studied the effect of the attractive range. The 
early studies [3, 15] considered perfect sticking of two monomers when their relative distance fell below the monomer 
diameter. In these works the authors did not use an interaction potential, but they examined the system frequently 
during its time evolution to identify agglomeration events (monomer-cluster or cluster-cluster) by calculating the 
relative distances of all monomers. After, e.g., two clusters had touched, the relative distances of all the monomers 
in the resulting cluster were "frozen", and the cluster was allowed to diffuse with a diffusion coefficient that had to 
be prescribed a priori. 

Herein, aggregate formation arises from monomer collisions that bind them through their interaction. We used two 
spherically symmetric intermonomer potentials: a model potential w™„ an d a potential u mm that arises from the 
integration of the intermolecular potential over the volume of two macroscopic bodies (e.g., two monomers). 

The model potential is short ranged, and it has a deep and narrow attractive well to model monomer binding 
without break-up. Furthermore, it tends smoothly to zero at r = r cut where r cut is a cut-off distance such that 
^cut ~ o o. This avoids the introduction of the so-called impulsive forces in the system [17]. It is attractive on a 
length scale much smaller than the monomer diameter cr. We chose the following analytical expression 
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where u a — u™^(er), r m ; n is the location of the potential minimum, and u m ; n = u mm ( r min)- The model potential 
depends on five parameters (r m ; n , cr, r cut , u a , and w m i n ) chosen in our numerical simulations as follows. The potential 
minimum is located at r m j n = 1.05er, where it evaluates to M m j n = — lOOfc^T. At monomer separation a the potential 
evaluates to u{a) = 60/csT with a steep gradient. For monomer separations in the range (0, cr] the potential is 
extrapolated linearly with the slope it has at r = a. Hence, for separations smaller than cr, the monomers feel a 
strong constant repulsive force equal to the force at r = cr. The cut-off distance is r cut = l.lcr. The model potential 
and its various parameters are shown in the inset of Fig. 1. 

The interaction potential may be obtained from the integration of the intermolecular interactions over the nanopar- 
ticle volumes as in Ref. [18]. We assume pairwise additivity of the intermolecular potential, continuous medium, and 
constant material properties. Elastic flattening of the monomer is neglected. The intermolecular potential is taken to 
be the Lennard- Jones potential 



ULj{r) = 4e 



where e is the depth of the attractive potential, the maximum attractive energy between two molecules, and cjlj the 
distance at which the potential evaluates to zero, the distance of closest approach of two molecules which collide with 
zero initial relative kinetic energy. The first term expresses (approximately) the repulsive part (Born repulsion), the 
second the attractive (van der Waals attraction). Integration of Eq. (11) over two equal-sized spheres of diameter cr 
yields, see, for example, Refs. [18, 19], an attractive part 
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to obtain 

u mm {r)=u^ m {r)+u<™{r) , (12c) 

where A = A-Kecr^jn 2 is the Hamaker constant and n the molecular number density in the solid. In the limit r ^> a 
the potential decays as u mm ~ — r~ 6 , i.e., the expected attractive van der Waals interaction energy between two non- 
charged macroscopic bodies is recovered. The Hamaker constant of a typical soot nanoparticle, for example n-hexanc, 
was estimated to be A — 2.38 • 10~ 19 J [20] with a corresponding cflj — 0.5949nm [21]. We used these values, along 
with a — 20nm, to render the interaction potential dimensionless. 

The repulsive part of u mm diverges for r — > er; it may, thus, cause numerical difficulties in the solution of the 
Langevin equations. We modified the potential at distances r ma t = 1.015cr (less than the position of the potential 
minimum at r min ~ 1.017<r) by extrapolating it linearly with the same slope it has at r mat until r = 0.995cr (where it 
evaluates to about 60/csT*). A similar matching condition was used to determine the coefficient of the model-potential 
linear term. We set u mm to zero at smaller monomer-monomer separations. This modification of the repulsive part 
is not expected to affect the dynamics of the system, since it involves monomer separations that are energetically 
unfavorable (and, hence, unlikely). Therefore, the intermonomer potential we use is 



(r) 



KZ( r ) =° if < r < 0.995<r , (13a) 
du 

(rmat - r) + u mm (r mat ) if 0.995ct < r < r mat , (13b) 

(r) = u mm (r) if r mat < r . (13c) 



dr 



We truncate the potential at distances r > 7<r where it is negligible with respect to the thermal energy, u mm (J o) I (ksT* ) 

io- 4 . 

Although neither potential diverges at separations r < er, at such distances monomers feel a very strong repulsive 
force; monomer separations below a are energetically unfavorable and their occurrence is extremely unlikely during 
the system dynamics. This justifies the identification of a with the hard-core monomer diameter. Moreover, since the 
two intermonomer potentials are much deeper than ksT* the sticking probability upon collision may be considered 
to be unity. 

In the following, we will use only dimensionless quantities; we will, thus, drop the tilde to simplify notation. 
Furthermore, unless specified otherwise, the results presented were obtained with w^j, Eqs. (13). Figure 1 presents 
and compares the two dimensionless radial potentials. 



III. NUMERICAL METHOD 



A. Numerical solver 



The numerical simulations were performed with the software package ESPRcsSo [22] , a versatile package for generic 
Molecular Dynamics simulations in condensed-matter physics. We used the Molecular Dynamics program with a 
Langevin thermostat. The Langevin thermostat was construed as formal method to perform Molecular Dynamics 
simulations in a constant temperature canonical ensemble, see, for example, Ref. [23]. It introduces a fictitious 
viscous force to model the coupling of the system to a thermal bath according to Eqs. (1) that ensures the system 
temperature fluctuates about the bath mean temperature. The molecular equations of motion with a Langevin 
thermostat are formally identical to the Langevin equations of interacting Brownian particles: the physical scales and 
their interpretation differ. 

The ESPResSo numerical solver uses the Verlet algorithm for the deterministic part, with a numerical error that 
scales at least like 0(St 3 ). The combination of the Verlet algorithm with the solver for the stochastic part (Langevin 
thermostat) yields an error estimate of 0(St 3 ) in the monomer positions and of order 0(St 2 ) in the velocities [24]. 
As a check of the numerical solver we simulated the motion of a single Brownian particle in a quiescent fluid. The 
simulations for the mean-square displacement (r\(t)), the mean-square velocity fluctuations (v 2 (t)}, and the velocity 
autocorrelation function (vi(t)vi(O)) agreed with the analytical expressions [14]. 

The initial state was created by randomly placing 71^ (0)V = 5000 monomers in a cubic box of size L with V = L 3 the 
box volume, and rioo (0) the initial monomer concentration. The box size was chosen to give the desired initial monomer 

1/3 

concentration according to L = (5000/noo) ' . We chose noo(0) = 0.01 corresponding to L ~ 80. The initial random 
placement of monomers could cause numerical problems if two or more monomers happen to be placed in almost 
overlapping positions, thus experiencing immediately a very strong repulsive force. We used a well-known technique 
of MD simulations to avoid these numerical instabilities by ramping up the repulsive force for r < 1 to its constant 
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final value during 800 time steps. The (dimensionless) simulation time step was chosen to be <5t s im = 1.25 x 10 -3 . 
After initialization, the system was evolved until a final time 3000, when the cluster concentration had decreased by 
almost two orders of magnitude. The results we show were obtained using the output of 10 simulations, each one 
with different initial monomer positions and zero initial monomer velocities. Finally, periodic boundary conditions 
were imposed. 



B. Cluster identification 



The identification of clusters is one of the most time-consuming tasks of post-processing the simulation results. 
The ESPResSo simulations return individual monomer positions and velocities. Unlike the previously mentioned 
works [3, 15], agglomeration events arc not identified during the time evolution of the system, but cluster formation 
is determined a posteriori. Sampling of the simulation output was performed at time intervals 5t sam = 2 (every 1,600 
simulation time steps). 

We resort to an approach based on graph theory. A cluster is a set of connected, bound, monomers. As both 
interaction potentials are very deep once two monomers collide and bind, they remain bound: no agglomerate break- 
up was noticed during our simulations. Hence, two monomers may be considered bound when their relative distance is 
less than a threshold distance c£ t hr- The choice of d t hr depends on the position of the minimum of monomer interaction 
potential r m ; n . For the simulations with the integrated Lennard- Jones potential we used c? t hr = 1-04, whereas for the 
model potential we used cZthr = 1.06. Small variations of c? t hr about these values did not affect the identification of 
clusters. Our definition of a cluster is reminiscent of the liquid cluster definition proposed by Stillinger and used in 
gas- liquid nucleation studies (see, for example, Ref. [25]). 

Computationally, the first step is the calculation of the distance matrix D between all monomers. For a box with 
periodic boundary conditions, the distance between two monomers is the distance between the i-th monomer and the 
nearest image of the j-th monomer [17]. For instance, the (ordered) distance between two monomers along coordinate 
I is 




= r{ - r} - L ■ nint -i— , l = x,y,z , (14) 



where "nint" is the nearest integer function. The periodicity of the box imposes a cut-off of L/2 on the maximum 
distance between two monomers along each axis. The three-dimensional distance 

a, = [E(^) 2 ] 1/2 ' (15) 

(=1 

is always a non-negative quantity. The distance matrix D is subsequently used to calculate the adjacency matrix A, 
defined as 

A^i 1 if A ^ rfthr ' (16) 
10 otherwise . 

The adjacency matrix is usually introduced in graph theory [12] as a convenient way to represent a graph uniquely 
Monomers in a cluster can be formally regarded as graph vertices (nodes), whereas the bonds due to the interaction 
potential become graph edges (links). The problem of cluster identification, given the distance matrix D, is then 
re-formulated as the identification of the connected components of a non-directed graph expressed by the (symmetric) 
adjacency matrix A. They can be determined using a standard breadth-first search algorithm [12]. Due to its speed 
and scalability, we resort to the implementation of its algorithm in the igraph library [26] . 



C. Cluster radius of gyration 



The radius of gyration is a geometric parameter used to characterize the size of fractal aggregates. It describes 
the spatial mass distribution about the aggregate center of mass. As such it is a static property that depends on the 
cluster mass distribution, and not on the diffusive properties of the cluster. For a cluster composed of k equal-mass 
monomers it becomes the root-mean-square distance of the monomers from the cluster center of mass [2] 

k 

R l = \il^-^CM? +a\ , (17a) 
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where the aggregate center of mass is 

k 

RcM^^r, , (17b) 

i=l 

and a\ a monomer characteristic size. Since we are interested in the power-law dependence of the radius of gyration 
on the number of monomers in a cluster even for small clusters we included a\ in the definition of the radius of 
gyration: otherwise Eq. (17a) evaluates to zero for a monomer. The additional term may be taken to be either the 
primary particle radius [2], a\ — cr/2, or the radius of gyration of a sphere [27], a\ — ^/3/5(er/2). The choice of ai 
as the monomer radius of gyration ensures that, if an average fractal exponent is used, the large k behaviour in the 
power-law scaling persists even for smaller clusters. Of course, the value of a\ (including a\ = 0) is irrelevant for large 
clusters. We used both choices for a\\ the results presented were obtained with a\ the monomer radius of gyration 
because this choice gave better agreement of the calculated and numerically estimated kernel elements, see Sec. V E. 
Nevertheless, the difference between the two choices was small. 

The calculation of R g according to Eq. (17a) must take into consideration the periodicity of the box. Since all 
monomer- monomer distances are known from Eq. (14), the position of monomers in the cluster with respect to a 
randomly selected monomer (say monomer 1) may be easily calculated; for instance, the j-th monomer position along 
coordinate I with respect to monomer 1 is 

r) = r[ + D\j . (18) 

Given the relative position of all monomers in the cluster, the position of the aggregate center of mass may be 
calculated via Eq. (17b), and the radius of gyration from Eq. (17a). We adopted a reference system centered on the 
randomly chosen monomer in the cluster, i.e., = (0,0,0). For this reference system rj = D[ .. This procedure is 
independent of the choice of the selected monomer, since the radius of gyration is independent of the cluster position 
in the box, nor does the procedure depend on the sign convention chosen for D[ . 



D. Collision kernel 



The package ESPResSo allows addressing each monomer individually at all times during the simulation, i.e. a 
permanent label (e.g., color) may be associated with each monomer. A monomer label allows the unequivocal identifi- 
cation of the cluster it belong to. Each cluster becomes an unordered collection of monomer labels where no monomer 
label is repeated: this amounts to the mathematical definition of a set. Viewing clusters as sets of monomer labels 
provides a computational method to investigate cluster collisions even for sampling times considerably longer than 
the simulation time step, as long as the aggregates do not break-up. To be specific, consider two clusters at time t. 
If during the (sampling) time interval (t, t + St sam ] they collide, they will be part of the same newly-formed cluster at 
t + <5i S am- An equivalent description of the collision is that the two monomer-label sets that identify the pre-collision 
clusters become proper subsets of the monomer-label set that identifies the cluster formed by their collision, and 
detectable at t + St sam . The number of collisions that occurred during (i, t + <5t sam ], and the clusters involved in the 
collisions, may be recorded simply by comparing different sets. 

The collision kernel Kij (rendered dimensionless by scaling it by <t 3 /ti) between an i-mer and j-mer during the 
sampling time interval is 

— = (2 - OijjKij— = (2 - d lJ )K lj n i (t)n J {t) , (19) 

where Nij{t) is the number of collisions between i- and j-mers that took place during the interval (t, t + St sam ], Ni(t) is 
the number of z-mers at time t, rii{t) = Ni(t)/V the cluster concentration, and Sij the Kroenecker delta. The collision 
kernel is estimated from Eq. (19) by treating it as the fitting parameter that minimizes (in a least-square sense) 
the distance between the number of detected collisions Nij(t)/(VSt sam ) and (2 — 5ij)Kijni(t)nj(t). 



IV. STATIC CLUSTER PROPERTIES 



Examples of 50-monomer clusters that survive at the end of one of the simulations are shown in Fig. 2. The top, left 
subfigurc is a snapshot of the system, whereas the other three are color-coded aggregates. The color code denotes the 
number of first neighbours of a monomer. Note the relatively compact and long, tubular structure of the generated 
aggregates. These two features, small-scale compactness and large-scale tubular structure, will be related to properties 
of the intermonomer potential in the following sections. 
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A. Cluster fractal dimension 

The fractal dimension of the clusters is determined from the statistical scaling law that governs the power-law 
dependence of the cluster radius gyration R g on cluster size k 

R g = ak^ d f , (20) 

where df is the average, time-independent fractal dimension of the aggregates, and a the fractal prefactor, occasionally 
referred to as lacunarity [28]. The fractal prefactor provides information on the packing of monomers [29]. The cluster 
radius of gyration was obtained by averaging over generated cluster configurations. Specifically, the average radius of 
gyration of a fc-monomcr aggregate is taken over all fc-monomer aggregates that had been recorded at least 200 times. 
This requirement aims at eliminating outliers. It is not particularly stringent since for each simulation 1500 system 
configurations (snapshots) were stored, and 10 simulations. Our results are not sensitive to reasonable choices of the 
threshold of minimum number of cluster occurrences: raising the occurrence threshold to 400 does not change the 
calculated fractal dimension. 

Figure 3 presents the calculated radius of gyration as a function of aggregate size. The double logarithmic plot 
shows that power-law scaling breaks down for cluster sizes k < 5. Hence, we set k = 5 the minimum cluster size 
in the fits. All cluster configurations with k > 5 were fitted to a single line leading to an average fractal exponent 
df = 1.62 ± 0.02, and a prefactor a = 0.242 ± 0.006. The prefactor re-expressed in terms of k = k g (2R g ) d f gives 
kg = 3.24. The calculated prefactor is on the high side of fractal prefactors reported in the literature for clusters 
generated by computer simulations [30], but closer to experimentally determined prefactors that indicate k g > 2 
[28, 30]. Inspection of the figure shows that clusters with k < 15 deviate significantly from the single-line fit. 
Simulation results are fitted better by considering two cluster population (and hence two different slopes) : wc identify 
small (k < 15) and large (k > 15) clusters. Refitting the data gives a fractal dimension d™ c = 2.25 ± 0.05 for small 
clusters and df — 1.56±0.02 for large clusters. The corresponding average prefactors, as well as the fractal dimensions, 
are reported in Table I. As before, the prefactor for the large clusters are on the high side of literature-reported values 
for computer-generated clusters, and closer to experimental values determined by angular light scattering, suggesting 
that the cluster generated in this study have a closely packed structure [29], see also Fig. 2. 

The presence of two cluster populations obeying scaling laws with different exponents may be attributed to different 
agglomeration mechanisms. Small clusters are generated mainly by monomer-cluster agglomeration; they are more 
compact and spherical than large clusters. This agglomeration process is similar to, but different from, diffusion-limited 
aggregation, as argued in Sec. IV B. Large clusters are mainly generated by cluster-cluster agglomeration. Therefore, 
it is expected, and numerically confirmed, that d™ c > df- Furthermore, the calculated fractal dimension of large 
clusters is comparable, but lower, to reported fractal dimensions of cluster-cluster agglomeration df ~ 1.7 — 1.8 [3, 31]. 
Further comments on the agglomeration process and the associated fractal dimensions are made in Sec. IV B. 

The existence of two cluster populations that arise from the predominance of different agglomeration mechanisms as 
agglomeration progresses suggests the calculation of a time-dependent average fractal dimension d f (t) . We calculated 
it by considering the instantaneous radius of gyration of fc-monomer clusters for each of the ten initial configurations: 
only clusters with k > 5 were included, and no occurrence threshold was imposed. As the number of clusters is 
considerably smaller than those used in the time-independent calculation we do not calculate first the instantaneous 
average radius of gyration and then fit the data (as done for the calculation presented in Fig. 3): instead wc fit all 
the data simultaneously (no averaging). The fractal dimension as a function of time is shown in Fig. 4 (top). The 
evolution of the fractal dimension may be linked to the kinetics of the small and large cluster populations. At early 
times, the system is almost entirely made up of small clusters (k < 15), Fig. 4 (bottom), hence df{t) ~ d™ c , whereas at 
late times the contribution of large clusters (k > 15) to the overall cluster population is dominant, hence df{t) —> d C f . 

These findings are in the spirit of the study by Kostoglou and Konstandopoulos [13] who used a distribution of 
fractal dimensions to characterize aggregate morphology. They showed that the mean fractal dimension relaxes to an 
a priori fixed asymptotic value specified by the dominant agglomeration mechanism. Our results on the importance 
of agglomeration kinetics on the evolution of df(t) are in agreement with earlier works, e.g., Refs. [3, 31] who suggest 
that the aggregate fractal dimension depends mainly on the dominant agglomeration mechanism. 

B. Cluster morphology 

As a measure of local compactness of the generated aggregates we calculated the cluster coordination number. 
The cluster coordination number, defined as the mean number of first neighbours of a monomer in a cluster, provides 
information on the openness of the aggregates and their compactness, and on the presence of cavities in their structure, 
i.e., their porosity. Hence, it is related to the fractal prefactor, Eq. (20); it characterizes the small-scale structure 
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of an aggregate, whereas the fractal dimension characterizes its large-scale structure. The coordination number 
varies between and 12 for spherical particles, reaching its highest value for hexagonal close packed (hep) or face 
centered cubic (fee) structures at a volume fraction of 0.74 [8]. The cluster coordination number is obtained during 
post-processing of the simulation results resorting again to the igraph library [26] . 

We calculated the mean cluster coordination number by averaging the coordination number of each cluster over all 
clusters at a given time. It is plotted in Fig. 5 as a function of time. At late times it reaches values higher than five, 
implying that the aggregates are relatively compact. 

The local, small-scale compactness of the aggregates, aggregates that are similar to those generated by the Langevin 
simulations of, e.g, Ref. [7], is related to the spherical symmetry of the monomer-monomer interaction potential. At 
early times when a monomer collides and binds to an aggregate it is free to slide, a motion induced by the thermal noise, 
and to reorient to maximize the number of contacts with other monomers. Since the potential is isotropic there is no 
angular restoring force to prevent sliding. However, short-time restructuring is hindered as the number of monomer- 
monomer contact points increase because the energy cost associated with stretching monomer-monomer bonds is high. 
Inspection of Fig. 2 shows that the minimum number of first neighrours in a stable cluster is three, suggesting that 
in three dimensions three contact points are sufficient to prevent monomer sliding. Hence, the aggregates generated 
herein are different from those generated by diffusion-limited aggregation where colliding monomers remain fixed at 
the initial point of contact. Becker and Briesen [10] showed that in the absence of a potential to prevent bending 
of monomer-monomer bonds the aggregate collapses to a more compact structure. The monomer reorientation at 
early times also suggests that the high fractal dimension d™ c we associated with monomer-cluster agglomeration 
mechanism arises from monomer rearrangement. In our simulations monomer reorientation occurred at very small 
time scales, smaller than the sampling time. After the short-time monomer reorientation clusters remain rigid till the 
next collision, as discussed Sec. IV C. 

An additional explanation for the compactness of aggregates generated by Langevin simulations, neglecting Brown- 
ian motion, was suggested by Tanaka and Araki [32] . They argued that in the absence of interparticle hydrodynamic 
interactions, in particular of squeezed-flow effects, the generated clusters tend to be more compact. 

Even though the clusters are locally compact, the large clusters generated at late times have a low fractal dimension 
dfj c < 2. Inspection of Fig. 2 (second and third clusters) suggests that the low, late-time fractal dimension is due to 
their tubular, elongated shape. The aggregates are not porous: they do not have holes nor cavities. The large-scale 
structure of the late-time (large) clusters is determined by the attractive range of the intermonomer potential. At 
late times when two (locally compact) clusters collide aggregate restructuring is limited to the monomers that are 
in contact: the attractive range of the potentials used in the simulation is too short to induce significant cluster 
re-organization, a process that would increase the fractal dimension. Videcoq et al. [7] showed that increasing the 
attractive range modifies the shape of the generated aggregates leading to more spherical aggregates. Hence, the 
late-time fractal dimension associated with cluster-cluster aggregation depends on the extend of restructuring due to 
the attractive potential range. This restructuring provides a possible explanation of the lower than expected fractal 
dimension for cluster-cluster aggregation. 

Additional simulations were performed to assess the sensitivity of the observed cluster morphology, and their rigidity, 
on the simulation time step. We simulated a 50-monomcr cluster (shown at the bottom right, Fig. 2) for 32000 time 
steps (tfinai = 40), with two different time steps, <5i s im and 5t s i m /4. We monitored the time evolution of the cluster 
radius of gyration R g (t), an average cluster property, and an instantaneous cluster property, the distance di ; i2(t) of 
two randomly chosen monomers (monomers 1 and 12). We found that both quantities fluctuated by approximately 
0.1%, independently of time step. These results confirm that cluster morphology does not depend on the choice of 
the time step, as long as it is chosen within reason. In fact, for considerably longer time steps the cluster breaks up. 
The fluctuations arise because monomer radial positions fluctuate, albeit slightly, about the potential minima. 

As a different measure of a possible dependence of cluster morphology on the simulation time step we compared the 
time-dependent adjacency matrices, Eq. (16), for these two simulations. The adjacency matrices were found to be all 
identical, another confirmation that cluster morphology is independent of the time step. Note that small fluctuations 
of monomer radial positions in the cluster are not reflected in the adjacency matrix (a binary matrix) since a distance 
threshold is used to convert the distance matrix, Eq. (15), to it. 

Thus, the morphology of the generated aggregates arises from a combination of kinetic effects, as induced by the 
thermal noise, and energetic effects, as determined by the attractive range of the interaction potential. A precise 
characterization of aggregate morphology requires both the fractal dimension and the coordination number (or the 
lacunarity) . 
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C. Cluster restructuring 

In addition to short-time monomer reorientation cluster restructuring may occur following a cluster-cluster collision. 
Cluster restructuring following contact of two clusters was investigated in a three-dimensional lattice model in Rcf. [33] . 
In our approach cluster restructuring does not have to be imposed as an additional feature of cluster dynamics, but 
it occurs naturally. The deep radial intermonomer potential locks distances between neighbouring monomers, but 
bonded monomers may slide over each other due to the thermal noise term. However, as mentioned earlier, monomer 
sliding is restricted by the number of monomer-monomer contact points. 

We use two statistical quantities as indicators of the modification of the structure of a cluster: the radius of gyration 
and the mean monomer-monomer distance. For a fc-mer the mean monomer-monomer distance D mm is the average 
monomer-monomer distance over all the monomers, 

1 k 

D m m = ^2 ^ij , (21) 

i,j 

where the monomer-monomer distance Dij was defined in Eq. (15). We monitored cluster restructuring by selecting 
a monomer at the beginning of the simulation and follow it as it collides with other monomers or clusters. In Fig. 6 
(left) we show D mm during a period when a collision occurs (t = 568) and the cluster size increases from k = 14 
to k = 43 . We notice that before and after the collision D mm is constant, a clear sign of cluster rigidity. When 
the collision between the two aggregates occurs, D mm increases, but on a time scale of a few monomer relaxation 
times it decreases, afterwards remaining relatively constant. We consider this behaviour a strong indication of cluster 
reorganization following the collision of two clusters. The cluster radius of gyration shows a perfectly analogue 
behavior, Fig. 6 (right). Note, however, that once cluster restructuring occurs after a cluster-cluster collision the 
resulting cluster retains its shape and it remains rigid until possibly the next collision. 

V. DYNAMIC CLUSTER PROPERTIES 

A. Cluster translational diffusion coefficient 

Equations (1) determine the dynamics of each monomer and indirectly the cluster diffusive properties. The cluster 
diffusion coefficient is calculated from the late-time dependence of the variance of the cluster center-of-mass position 
as a function of time [14] 

(SR 2 CM (t)) = ([R C mW - (RcmW)] 2 ) ^^6D k t . (22) 

A different set of simulations was performed to determine the diffusion coefficient of a few selected clusters. As noted 
in Sec. IV C in the absence of collisions clusters are rigid: the dynamical properties presented in this section refer 
to perfectly rigid clusters. A selected cluster was placed in the simulation box with its centre-of-mass position at 
Rcm = (0,0,0) and with zero monomer velocities. The box size was chosen to be L = 10,000, a size much larger 
than the box size used for the agglomeration simulations. A large box is necessary because the aggregate centre of 
mass has to be tracked for a (potentially) long time, and, due to the periodic boundary conditions, no displacement 
of the aggregate from its initial position larger than L/2 is admissible. For these simulations the cluster never crossed 
the boundary of the box. 

The center of mass of a single cluster was tracked up to t = 400 for clusters with k = 4, 10, 50, 98. Averages were 
performed over 800 trajectories, each trajectory starting with an identical cluster: different cluster trajectories arise 
from different realizations of the stochastic noise. Figure 7 (top) shows (6RQ M (t)) for a 50-monomer cluster, chosen 
to be the cluster at the bottom right of Fig. 2. The time dependence of the ensemble-averaged, mean-square cluster 
displacement quickly becomes linear, in agreement with Eq. (22). The inset of Fig. 7 (top) magnifies (on a double 
logarithmic scale) the early-time behavior of (SR^, M (t)). For t < 1 cluster motion is ballistic, and the mean-square 
displacement exhibits a power-law dependence on time, (6R^ M (t)) oc V with 7 = 3. This is the expected behaviour 
for a single Brownian monomer with zero initial velocity [14]. Similar early-time ballistic motion was found for a 
cluster with k = 4. 

The results of the numerical simulations for the four clusters are summarized in Table II, where we also report 
the cluster radius of gyration (obtained from Fig. 3). They may be interpreted by considering theoretical estimates 
of the diffusion coefficient of a fc-mer. The Stokcs-Einstcin diffusion coefficient of a single spherical monomer, in 
dimensionless form, scaled by D* = ksT* /(/3imi), is [14] 



Di=T 



(23) 
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The diffusion coefficient of a fc-mer of mass M k = km\ may be expressed as a generalization of Eq. (23). Let f3 k be 
the average friction coefficient of a monomer in the fc-clustcr. Then, the cluster drag term in Eq. (1) may be written 
as (3 k M k = fc/3fcmi, and the cluster diffusion coefficient generalizes to 

Dk Dl kf3 k ~ Dl kr) k ' ^ 

Equation (24) manifestly shows the importance of the ratio of the average friction coefficients rj k = (3 k , rj k being 
the average shielding factor of a monomer in a fc-cluster. In fact, Filippov et al. [2], in a different context, used the 
cluster shielding factor to describe how monomer shielding affects heat transfer to an aggregate. Herein, the cluster 
shielding factor is connected to cluster mobility. 

The numerically determined diffusion coefficients behave as D k oc 1/fc with k spanning almost two orders of 
magnitude. The estimated value of (3 k equals the monomer friction coefficient /?i within a few percents, Table II 
(/?i = 1 in our units). The small deviations are attributable to the limited number of simulated stochastic trajectories. 
Our Langevin simulations suggest that the cluster diffusion coefficient is inversely proportional to the cluster mass and 
to the monomer friction coefficient. This behaviour is a direct consequence of neglecting shielding of inner monomer 
by outer monomers in an aggregate. The cluster accessible area becomes the sum of the monomer accessible areas, 
since for r] k = 1 the total friction coefficient of a fc-mer is 3ir[if<jk. For this reason we refer to these clusters as ideal 
clusters, with respect to their dynamical properties. For ideal clusters the average monomer shielding factor is unity 
irrespective of the state of aggregation of the monomer. This is an inherent property of all simulations of monomer 
agglomeration via the unshielded monomer Langevin equations Eqs. (1). 

The cluster diffusion coefficient (in the continuum regime) is frequently expressed in terms of a mobility radius R m , 
the radius of a sphere with the same friction coefficient as the aggregate [34] . In the dimensionless units used in this 
work, the mobility radius determines the cluster diffusion coefficient by 

D k = D!^- . (25) 

Comparison of Eqs. (24) and (25) gives an expression for the mobility radius in terms of the (average) cluster (per- 
unit-mass) friction coefficients f3 k and or the shielding factor rj kl 

k j3 k k 

Rm= 2j 1 = 2 1]k ■ (26) 

For the ideal clusters generated in this work R m = k/2. Table II presents the mobility radius calculated either from 
the numerically determined diffusion coefficient [Eq. (25), fourth row], or by evaluating it directly from Eq. (26) with 
f] k = 1 (fifth row). A comparison of these two values provides an indication of the numerical error of our simulations: 
the percentage shown in the fifth row (in parenthesis) quantifies the difference. 

The cluster diffusion coefficient is frequently estimated by taking the aggregate mobility radius R m approximately 
equal to the cluster radius of gyration R g . Our simulations show that the radius of gyration is of the same order of 
magnitude for cluster sizes in the range k = 4 to k = 98, whereas the mobility radius is considerably higher. Since 
the radius of gyration severely underestimates the mobility radius, the diffusion coefficient of ideal clusters is lower 
than the diffusion coefficient of clusters with R m ~ R g . In ideal clusters the average monomer shielding factor is 
unity, whereas monomers in non-ideal clusters are at least partially shielded, resulting in higher cluster mobility, cf. 
Eq. (24). 



B. Cluster thermalization 



For a single Brownian monomer, the late-time mean-square velocity fluctuations (dimensionless, scaled by fc^T* /mi) 
tend to [14] 

(Svl 00 )= lim([v 1 (i)-(v 1 (i))] 2 )=T , (27) 

t— too 

where Vi(t) is the monomer instantaneous velocity. For a fc-mer the previous expression, a manifestation of energy 
equipartition and a consequence of the FDT, may be generalized to 

(SVCM IOO ) = ([VCMM - (VcMit))} 2 ) = <*«l,oo> l > ( 28 ) 
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where the cluster centrc-of-mass velocity is Vcm(£) — X)i v i/^j with Vj the z-th monomer velocity. As for a Brownian 
monomer Eq. (28) expresses energy equipartition for a fc-monomer cluster. In Fig. 7 (bottom) we plot {5V^, M {t)) as 
a function of time for the 50- monomer cluster used in the previous section. As shown, Fig. 7 (bottom), (SV^, Moo (t)) 
tends at late times to 0.03 (up to noise fluctuations), i.e., the theoretical value for ideal clusters, Eq. (28), for the 
simulation parameters T = 0.5 and k = 50. 

The main result of Sees. VA and VB is that a monomer in an ideal fc-cluster generated according to Eqs. (1) and 
(2), i.e., neglecting monomer shielding, feels the same friction coefficient as an isolated monomer (/3& = /?i). An ideal 
fc-monomer aggregate has the same diffusive properties as a monomer sphere of k times the mass of a single monomer 
(Dk = Di/k). The cluster diffusion coefficient differs from the monomer diffusion coefficient only due to the larger 
cluster mass, and not, in addition, due to the decrease of the average monomer shielding factor. 



C. Cluster rotation 



We investigated whether ideal clusters have a preferential orientation by examining the distributions of their average 
Euler angles. As argued, clusters in the absence of collisions may be treated as rigid bodies. Rigid body rotation 
in three dimensions may be described by the three Euler angles 6, <j> and ip [35]. We evaluated them during post- 
processing by eliminating cluster translational motion via rigidly translating the cluster centre of mass to the origin 
of the computational-box coordinate system at all times. The Euler angles may, then, be calculated by recording 
the time-dependent positions r^(i), rs(i) and rc(t) of three non-coplanar monomers. We define a 3 by 3 matrix 
X(0) = [i"a(0), i\b(0), rc(0)] containing the initial positions of the three reference monomers and a matrix X'(i) = 
[r^(i), r B (i), rc(i)] with their positions at time t. The rotation matrix A such that X' = AX is 

A(t) = X'(t)X T (0) [x(0)X T (0)l ~* , (29) 

where the superscript T denotes matrix transpose. Once the rotation matrix A is known, the Euler angles may be 
determined. Since they are not uniquely defined, their values depending on the order of the three rotations, we chose to 
calculate them via the algorithm presented in Rcf. [36]. The range of the Euler angles so determined is: <f) an d i> range 
in the interval [— tt,tt], whereas 8 lies in the interval [— tt/2, 7r/2]. For uniform random rotation matrices [37], both ip 
and <f> are random variables with a uniform probability distribution in the interval [— 7r,7r]: hence (cf>) = (ip) = and 
(Stf it)) 1 / 2 = ([<f>(t) - (0(t))] 2 ) 1 /2 = (^2(^1/2 = n ^ _ 1 gl on the other hand, the Euler angle 6 is distributed 
according to [37] 

6 = arccos(l - 2Z(0, 1)) - | , (30) 

where Z(§, 1) is a random variable with a uniform probability distribution in the interval [0, 1]. The corresponding 
averages, obtained numerically, are (0(f)) = and ((50 2 (t)) 1 ^ 2 ~ 0.68. 

Figure 8 (top) shows the mean Euler angles calculated for an ensemble of 800 identical clusters containing 10 
monomers each. Their mean values fluctuate about zero. On the bottom diagram we show their mean-square 
fluctuations: they rapidly tend to the expected values for random rotation matrices. Another calculation of (89 2 ) for 
an ensemble of 800 identical cluster with k — 50 shows saturation to the same value but on a longer time scale. The 
rotational thermalization time of an ideal fc-cluster depends on its total mass. Hence, at late times an ideal cluster 
does not have a preferential orientation. 



D. Time dependence of cluster number: The agglomeration equation 



The instantaneous total number of clusters, N^t) = Nk(t), quantifies the progress of agglomeration. It may 
be calculated via the numerical solution of the agglomeration equation [34] 

< ^ = \ ^2 K a n i n j ~ n k^K ik m . (31) 

i+j — k i 

For definitions and scalings see Section HID. In this section the collision kernel Kij is scaled by K* = D*a = /3ict 3 , 
see Eq. (6). 

The standard treatment of the collision kernel for fractal aggregates in the continuum regime gives [34] 

= 4tt( Di + Dj){Ri +Rj) , (32) 
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where we dropped the subscript g in the radius of gyration of an i-mer, Ri = R g ^, to simplify notation. Equation (32) 
is often expressed in terms of the aggregate volume of solids Vi oc (Ri) dj [38] since the volume of solids is conserved 
during agglomeration. Moreover, the diffusion coefficients in Eq. (32) are given by Eq. (25) with the mobility radius 
being the radius of gyration. The kernel then reads 

*r + u'* + »r) , (33) 



where the superscript Sm refers to the so-called Smoluchowski kernel. The Smoluchowski kernel is not expected to 
model our numerical results accurately because it is based on diffusion coefficients given by Eq. (25) with R m = R g , 
an equality not respected by our simulations. Instead, the numerical simulations may be used to derive a modified 
kernel appropriate for the reproduction of the numerical results. Such a kernel is obtained by using the numerically 
determined cluster diffusion coefficients, Eq. (24) with = /?i in Eq. (32) For the radius of gyration, we use the 
fitted value for k > 5 (and the data reported in Table I) to obtain 



Ri 



simulation data, i < 5 , 

a mc i 1/d 7\ 5<z<15 , (34) 

a cc i 1/d T, 15 <i . 



Furthermore, non-continuum effects arising from monomer-monomer collisions, and dependent on the monomer mean 
free path, have been shown to be important in combustion-generated nanoparticlc agglomeration [38]. They may be 
accounted for by introducing the (dimensionless) Fuchs correction factor (3p in the kernel: for the explicit expression 
of /3f consult Rcf. [38]. Hence, the kernel appropriate for the Langevin simulations (hereafter referred to as Langevin- 
Dynamics kernel) reads 

Kf/ D =^D 1 + -^j (Ri + Rj) p F . (35) 

The late-time dependence of the total cluster number may be expressed in terms of the kernel homogeneity exponent 
A [39], whereby a kernel is a homogeneous function of order A if K^^j = r y x Kij. Then, the asymptotic time-decay 
of the total cluster number is iVoo(£) ~ t -1 ^ 1- ^. For the Smoluchowski kernel A — and, therefore, ~ t~ x , 
whereas for the Langevin-Dynamics kernel A = (l/dj ) — 1, leading to <~ i~ - 74 . 

The results of the simulations (with both potentials) for the total number of clusters and the numerical solution of the 
agglomeration equation with K Sm and K LD are shown in Fig. 9. The numerical solution of the agglomeration equation 
with the Langevin-Dynamics kernel (short dashes) shows good agreement with the simulations (dots, simulated 
intermonomer potential) The early-time agreement is attributed to the Fuchs correction factor, an observation also 
made in Ref. [38]. 

The time decay of the total cluster number was fitted to a power law <~ at times 2500 < t < 3000. For the 
Langevin-Dynamics kernel we find £ = 0.77, a value close to the exponent determined from the numerical simulations 
(£ = 0.78 and £ = 0.79 for model and the simulated potentials, respectively) and to £ = 0.74 expected for homogeneous 
kernels. The Smoluchowski kernel exponent is considerably different, £ = 1. The slight difference between calculated 
and theoretical exponents is attributed to the slow approach to the asymptotic limit of the agglomeration equation: 
the asymptotic limit is reached at times about one order of magnitude longer than the duration of the simulations, 
primarily due to the Fuchs factor. Finally, we note that, as expected, the Van der Waals attraction enhances the 
agglomeration rate (compare the simulation results with the two potentials), as noted in Rcf. [7]. 



E. Collision kernel elements 



Collision kernel elements may be obtained by ensemble averaging Eq. (19) over the 10 initial conditions, i.e., 

= {2-5 ij )K ij {n i {t)n j {t)) . (36) 



(Nij(t)) 



We found that the average of the nonlinear term decouples, (riirij) = (ni)(rij), to the accuracy of our simulations. 
We recorded cluster collisions every St sam , to avoid counting ternary collisions. We had enough data to calculate 
only a few kernel elements K t j for low i and j indexes. Figure 10 shows an example of the fitting procedure to 
determine K13 from (N13) / (VSt sanl ) and (ni)(n 3 ). The calculated kernel elements are reported in Table III, and 
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they are compared to the analytical values of . Kernel elements are scaled to the diagonal of the Smoluchowski 
kernel Kf™ — 8ttDi. We note an excellent agreement for the {(1, 1), (1, 2), (2, 2)} elements. These kernel elements are 
important in the early-time behaviour of Noo(t), and therefore they are responsible for the agreement of the numerical 
simulations with the numerical solution of the agglomeration equation at early times, as shown in Fig. 9. 



VI. SUMMARY AND CONCLUSIONS 

The static structure and diffusive properties of aggregates formed in a quiescent fluid by the collision, and subsequent 
binding, of spherical monomers were investigated via Langevin dynamics, in the limit of small Knudsen number 
(continuum regime). The Langevin equations of motion of a collection of diffusing and interacting monomers were 
solved numerically using a package for generic molecular dynamics simulations with a Langevin thermostat. Two 
intermonomer interaction potentials were used: a model potential and a potential that arises from the integration 
of the Lennard- Jones intermolecular potential over the volume of two monomers. Both potentials arc spherically 
symmetric and rapidly decaying. Aggregates were identified during post-processing by considering them to be the 
connected components of a non-directed graph. 

The static structure of the generated aggregates was described in terms of their average fractal dimension and 
cluster coordination number. We found that the aggregate fractal dimension varied with time from an early-time value 
characteristic of monomer-cluster agglomeration, (df — 2.25 ± 0.05), as determined by local monomer rearrangement 
and restructuring, to a late-time value characteristic of cluster-cluster agglomeration (df = 1.56 ± 0.02), a value 
dependent on the attractive range of the intermonomer potential. The time dependence of the fractal dimension was 
linked to the dynamics of two cluster populations, small clusters (k < 15) at early times and large clusters (k > 15) 
at late times. The average cluster fractal dimension, thus, was related to the dominant agglomeration mechanism. 
The generated aggregates had a cluster coordination number, defined as the mean number of first neighbours of a 
monomer in a cluster, of more than five (at late times) suggesting that the clusters were compact and not porous. The 
generated clusters were long, compact, and tubular with a high mean coordination number and low fractal dimension. 

We argued that these two salient features of aggregate morphology, small-scale local compactness (evident at 
the early stages of the agglomeration process) and longer-scale tubular structure (evident at later stages of the 
agglomeration process), were a consequence of properties of the monomer-monomer interaction potential. The small- 
scale structure is determined by the isotropy of the potential that allows short-time local reorientation of bonded 
monomers induced by the thermal noise. Monomers are free to rearrange to maximize their interaction with other 
monomers, since there is no angular potential to hinder bending of monomer-monomer bonds, under the constraint 
that monomer-monomer bonds are not stretched. The large-scale structure is determined by the potential interaction 
range. Colliding, locally compact, clusters rearrange locally at the point of contact: the attractive interaction range 
of the potential used in the simulations was too short to induce larger-scale rearrangement. Increasing the attractive 
range would lead to more spherical aggregates. After the initial rearrangements the aggregates remained rigid till the 
next collision. 

The dynamic cluster properties were analyzed in terms of the cluster translational diffusion coefficient. We found 
that the diffusion coefficient of a fc-mer scaled with cluster size as Dk oc fc _1 : aggregates diffuse like massive monomers. 
Furthermore, the average (per-unit-mass) friction coefficient of a monomer in a fc-monomer cluster was found to equal 
the (per-unit-mass) friction coefficient of an isolated monomer. Hence, the friction coefficient of a monomer in a 
cluster was determined to be independent of its state of aggregation, an approximation referred to as free draining 
approximation: the average monomer shielding factor of the generated clusters was unity. We argued that this diffusive 
behaviour is a consequence of the absence of shielding in the Langevin equations. This is an inevitable consequence of 
the usual application of Langevin simulations, unless a (time dependent) shielding coefficient is explicitly introduced 
in random force term of the Langevin equations. The shielding coefficient would, as a consequence of the Fluctuation 
Dissipation Theorem, modify the monomer friction coefficient. We, thus, referred to these clusters as ideal clusters, 
with respect to their transport properties. Similarly, the generated clusters did not have, on average, a preferred 
orientation. 

We also calculated numerically and compared to analytical predictions kernel elements Kij for low indexes. An 
extensive numerical investigation of the kernel elements is beyond the purpose of the present study, but it can be 
performed using the methodology described herein. 

We conclude that aggregates generated by unshielded Langevin equations of motion of monomers interacting via a 
central, rapidly decaying potential are on small scales, locally compact, on larger scales tubular, and they diffuse as 
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massive monomers. 
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TABLE I. Overview of cluster static properties. 





Fractal dimension 


Prefactor (R g = ak 1/d f) 


Prefactor (k = k 3 (2R g ) d f) 


All clusters 
k < 15 
15 < k 


d f = 1.62 ±0.02 
dj c = 2.25 ± 0.05 
df = 1.56 ±0.02 


= 0.242 ±0.006 
a mc = 0.386 ± 0.009 
a cc = 0.218 ±0.007 


k g = 3.24 
fc™ c = 1.75 
kf = 3.65 



TABLE II. Overview of cluster diffusive properties. The percentage in parenthesis (fifth column) is the relative difference of 
the two mobility radii reported in the fourth and fifth columns. 







k = 4 


fe = 10 


fe = 50 


fe = 98 


D k 


Simulations, Eq. (22) 


1.30- lO" 1 


4.92 • 10^ 


9.48 • 10~ d 


4.84 ■ 10~ d 


Rg 


Simulations, Fig. (3) 


7.5- 10 _i 


1.14 


2.12 


4.70 


Pk 


Simulations, Eq. (24) 


9.6- lO" 1 


1.01 


1.05 


1.05 


Rm 


Simulations, Eq. (25) 


1.92 


5.08 


26.35 


51.56 


Rra 


Eq. (26) with t] k = 1 


2 (4%) 


5 (1.6%) 


25 (5.4%) 


49 (5.2%) 



TABLE III. Comparison of kernel elements (expressed in units of 8ttDi). 



(i,3) 


K$ u [Eq. (35)] 


Kij [Eq. (36)] 


(1,1) 


0.213 


0.214 


(1,2) 


0.294 


0.298 


(1,3) 


0.309 


0.397 


(1,4) 


0.309 


0.449 


(2,2) 


0.317 


0.317 


(2,3) 


0.303 


0.349 




FIG. 1. Model and simulated monomer-monomer interaction potentials. 
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1 2 5 10 20 50 100 200 500 

k 

FIG. 3. Average radius of gyration as a function of cluster size. Linear fits performed on a double-logarithmic scale for k > 5 
(long-dashed line), 5 < k < 15 (short-dashed line) and k > 15 (solid line) (see, also, Table I). 
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FIG. 4. Top: Aggregate fractal dimension as a function of time. Bottom: Total number of clusters (solid), number of small 
clusters k < 15 (short dashed), and number of large clusters k > 15 (long dashed). 
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FIG. 5. Mean cluster coordination number. 
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FIG. 6. Top: Mean monomer-monomer distance of a selected cluster. Bottom: Radius of gyration of the same cluster. 
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FIG. 7. Top: Mean-square cluster displacement averaged over 800 trajectories (crosses) and linear fit (solid line) of a 50- 
monomer cluster (shown in Fig. 2, bottom right). Inset: Early-time behaviour of (SRom) (crosses), linear fit (solid line) and 
power-law fit, 7 = 3 (dashed line). Bottom: Velocity fluctuations {SVcm) (crosses) and analytical expression, Eq. (28) (solid 
line). 
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FIG. 8. Top: Mean Euler angles (6) (solid line), {(j>) (dashed line) and (if)) (dot-dashed line), averaged over 800 trajectories, for a 
10-monomer cluster. Bottom: Standard deviations (<5# 2 ) 1//2 (solid line), (5(f) 2 ) 1 '' 2 (long-dashed line), and (Sip 2 ) 1 ^ 2 (short-dashed 
line) for the same 10-monomer cluster, and (SO 2 ) 1 ^ 2 (dot-dashed line) for a 50-monomer cluster. The horizontal lines are the 
standard deviations of random rotation matrices. 
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FIG. 9. Total number of clusters as a function of time calculated with the simulated intermonomer potential (diamonds), model 
potential (filled circles), numerical solution of the agglomeration equation with the Smoluchowski kernel (long-dashed line) and 
with Langevin-Dynamics kernel (short-dashed line). 
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FIG. 10. Determination of kernel element K13: {Nis)/St diamonds, 2Ki^{ni){nz) with the fitted value of K13 circles. 



